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ABSTRACT 

We simulate the evolution of halo wide binaries in the presence of th e MAssive Compact Hal o Objects (MA- 
CHOs) and compare our results to the sample of wide binaries of Chan ame & Gouldl (120031) . The observed 
distribution is well fit by a single power law for angular separations, 3 "5 < Af? < 900", whereas the simulated 
distributions show a break in the power law whose location depends on the MACHO mass and density. This 
allows us to place upper limits on the MACHO density as a function of their assumed mass. At the 95% confi- 
dence level, we exclude MACHOs with masses M > 43 M Q at the standard local halo density pu- This all but 
removes the last permitted window for a full MACHO halo for masses M > 10~ 7 5 M Q . 

Subject headings: dark matter — Galaxy: halo — methods: numerical — stars: binaries — stars: kinematics 
— stars: statistics 



1. INTRODUCTION 

After formation, wide binaries retain their original orbital 
parameters except in so far as they are affected by gravi- 
tational encounters, or either one or both binary members 
evolve off the main sequence. These systems can therefore be 
used to probe inhomogeneities of Galactic potential that may 
be due to black holes, low luminosity stars, molecular clouds, 
or other objects, because their low binding energies are easily 
o vercome by gravitational p e rturb ations (Heggi eil975l) . 

Bahcall. Hut & Tremaine (1985) were the first to apply 
this principle, using wide binaries in the Galactic disk 
to inv estigate disk dark matter in the Sola r neighbor- 
hood. IWeinberg. Shapiro & Wassermanl (119871) refined this 
approach by incorporating several effects that were previously 
ignored. In particular, they showed that gravitational encoun- 
ters do not in general induce a sharp cutoff in the binary 
semi-major axis distribution, but rather generate a break in its 
power law. Although there have been extensive investigations 
of disk binary systems, they have not yielded strong conclu- 
sions about the disk potential. This is partly due to the rela- 
tively small size of the disk binary samples available, partly 
to the fact that they do not have good sensitivity for separa- 
tions > 0.1 pc, and partly because of the intrinsic complexity 
of the disk potential. In principle, halo binaries could be used 
to search for dark matter, but this has never been done, pri- 
marily because there were no halo-binary samples adequate 
to the task. 

However, halo dark matter has been in vestigated by sev- 
eral other techniques. lLacev & Ostrikerl (Q985) suggested 
a mechanism for disk heating by supermassive black holes 
and discussed that black holes with M > 10 6 M Q could 
destroy the disk. Null results of a search for "echoes" 
of gamma-ray bursts induced by gravitational lensing con- 
strained halo dark matter in the m ass rang e M ~ 10 6 5 - 
10 85 M Q JNemiroffetalJ 119931) . iMoord (119931) argued 
that massive black holes M > 10 3 M Q could disrupt low- 
mass globular clusters, which would imply an upper limit 
M < 10 3 M Q . However, this argument is somewhat sen- 
sitive to assumptions about the initial population of glob- 
ular clusters. Finally, microlensing exp eriments by the 
MACHO collaboration (Alco cket alJl200l . the EROS col- 



laboration iAfonso >s\ alJ|2Q 03l). and t he two collaborations 
working together (Alcocketal. 1998) found that MAssive 
Compact Halo Objects (MACHOs) with 10" 7 5 M Q < M < 
30 M<^ cannot account for the mass of the dark halo. 
iDe Ruiula. Jetzer & Massol l[lj?92) argued that baryonic dark 
matter with M < 10~ 7 M Q would have evaporated away in a 
Galactic time scale. 

The publication of the \i > 0."18 yr _I prop er motion lim- 
ited New Luyten Two Te nths (NLTT) catalog (L uvtenlll979t 
iLuvten & Hughes! 119801) has vastly increased the pool of 
available data on binary systems, but the short color baseline 
of the photographic photometry in NLTT rendered halo/disk 
discrim inati on extremely difficult. However, iGould & SalirrJ 
J200l and ISalim & Gouldl (|2003) revised NLTT (rNLTT) 
with improved astrometry and photometry for the 44% of the 
sky covered by the intersection of the Second Incremental Re- 
lease of the Two Micron All Sky Survey and the first Palomar 
Observatory Sky Survey. Chaname & Gould (2003, hereafter 
CG) then constructed a homogeneous catalog of binaries from 
rNLTT and classified each entry as either disk or halo. 

In this paper, we use Monte Carlo simulations to evaluate 
the effect of MACHOs on the halo binary distribution and 
compare these predictions with the observed halo binary sam- 
ple of CG by means of a likelihood analysis. We briefly de- 
scribe our CG halo binary sample in §|2]and make simple ana- 
lytic estimates of the effects of perturbers in §|3] Detailed jus- 
tifications of our assumptions and our Monte Carlo algorithm 
are presented in § 0] We present the results of our numeri- 
cal simulations in §|5]and the likelihoods of halo dark-matter 
models in §[6] Finally, we summarize our results in §0 

2. THE CG HALO WIDE BINARIES 

Figure^]shows the halo binaries from the CG catalog upon 
which the current work is based. This binary distribution 
function is well approximated by a single power law from 
AO = 3 "5 to the catalog limit at A9 = 900". However, since 
the lower threshold of completeness is not precisely estab- 
lished, we consider several different lower limits. There are 
respectively (90, 68, 59, 47) binaries in the samples with 
lower limits (3"5, 5"5, 7"5, 10"5). For AO < 3"5, the sample 
is not complete, and the deviation from a power-law distribu- 
tion is clearly shown as triangles in Figure^ However, since 
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FIG. 1. — Halo bi nary distribution function from the catalog of 
Chaname & Gould 12003). Samples of wide binary systems are complete 
up to 900" while the lower limit in angular separation is not precisely estab- 
lished. The circles represent actual counts in six equal logarithmic bins over 
3"5 < A8 < 900". The triangles show three \" bins for A9 < 3"5 and 
are rescaled to account for the smaller bin sizes relative to the circles. Vari- 
ous lines represent fits for subsamples with different lower limits in angular 
separation. Error bars represent one standard deviation cr = ra"'/ 2 /ln 10. 



the fit is insensitive to the choice of lower limit (A6> > 3. "5), 
we adopt a lower threshold of A9 = 3 "5. Since severe incom- 
pleteness clearly sets in at AO < 3 "5, one might be concerned 
about more modest incompleteness just above our adopted 
boundary AO = 3 "5. However, as we discuss in §[6] incom- 
pleteness near AO ~ 3. "5 would act to lessen the likelihood 
difference between a halo dark-matter model and the obser- 
vations. Therefore, our procedure is conservative. 

Assuming that the initial binary distribution is characterized 
by a power law, the slope of which might be different from the 
final observed one, the observed binary distribution can yield 
significant constraints on halo dark matter. Since binary sys- 
tems with large semi-major axes, a > 0.1 pc, are more liable 
to be disturbed by perturbers than those with small a, the de- 
viation (or lack thereof) from a power-law distribution at wide 
angular separations provides a test of halo dark-matter mod- 
els. 

3. THEORETICAL EXPECTATIONS 

We begin by considering two extreme regimes; the tidal 
regime (b m i n ^> a) and the Coulomb regime (b^ <C a), where 
a is the binary semi-major axis and, 
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is the typical minimum impact parameter for perturbers of 
mass M, density p, and velocity v over the duration T that the 
binary is subjected to perturbations. We evaluate equation Q 



using the standard local halo density pn = 0.009 M Q pc 3 
(Bahc all & Soneiral H980). We adopt T = 10 Gyr. Note that 
while halo stars (and so halo binaries) are several Gyr older 
than this, T represents the time since the binaries were as- 
sembled into the Galactic halo and not the time since their 
formation. 

In the tidal regime, the perturbations are dominated by the 
single closest encounter, which yields a change of relative ve- 
locity, 

Av~-2-a. (2) 
b L v 

By evaluating Av at b = b m m and equating the result with the 
internal velocity of the binary system, v 2 = Gm/a, we obtain 
an estimate of transition separation, 
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where we adopt m = 1 M as the typical total mass of the 
binary system. The corresponding transition orbital period is 
P t ~2Myx(p/p H )-\ 

Binary systems whose semi-major axes are less than a, or 
equivalently whose orbital periods are shorter than P, are in- 
sensitive to perturbations. Note that the transition separation 
in the tidal regime is independent of the mass and veloc- 
ity of perturbers. This regime approximately applies when 
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In contrast to the tidal regime, which is dominated by the 
single closest encounter, the perturbations in the Coulomb 
regime (named for the Coulomb logarithm, In A = ln(a/ b m i n )), 
are described by continuous weak gravitational encounters, 
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where the factor 2 in front accounts for the independent per- 
turbations of each component of the binary. The transition 
separation is then, 
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and so is roughly inversely proportional to the perturber mass. 

4. NUMERICAL SIMULATION 

In this section, we summarize our Monte Carlo algorithm, 
which evaluates the effect of all perturbations using a simple 
impulse approximation and ignores such effects as large-scale 
tides and molecular clouds. We include the con tribu tion from 
the ionized binaries, although as we discuss in § !4.4l and show 
more fully in the Appendix, disrupted binaries diffuse to sep- 
arations well beyond the observational range of interests in a 
time very short compared to T. 

Since previous work, nota bly that of 
Weinberg. Shamro & Wasserman (119871) . has focused 
considerable effort on incorporating the above-mentioned 
effects, we first justify our decision to ignore them. It is 
primarily the fact that we are considering halo binaries 
whereas previous workers were investigating disk binaries 
that accounts for the difference in importance of these effects. 



3 



4.1. Impulse Approximation 

In the Coulomb regime, the biggest relevant impact param- 
eter is of order a. Since the perturbers are moving much faster 
than the binary components, the impulse approximation well 
describes encounters. 

In the tidal regime, the perturbations are dominated by the 
closest single encounter, whose crossing time is, 
2h ■ 

*c — 

Since, for masses M < 10 8 Mq, the crossing time is signifi- 
cantly less than the transition orbital period P, ~ 2 Myr, the 
impulse approximation is always justified. 

4.2. Disk and Galactic Tides 

When halo binaries pass through the plane of Galactic disk, 
differential gravitational attraction s give rise to tidal effects, 
also known as disk shocking (Binnev & Tremaine 1987). 
Since the mean velo city at which ha lo binaries cross the disk 
is v, ~ 100 km s -1 (Popowski & Gould 1998a b), the change 
in the z-component of velocity of halo binaries is Av z ~ 
AnGY^a / v, where £ = 40 M Q pc" 2 is the su rface density of the 
Galactic disk in the Solar neighborhood (Zheng et al. 2001, 
and references therein). Equating this to the internal velocity 
of the binary system, v 2 = Gm/a, yields, 

/ 2 \ V3 

/ mvf \ 

flcrit = I 167r 2 GS 2 ) - 2 P C , (Disk Tides). (8) 

For the most favorable (i.e., prograde) orbits, the binary sys- 
tem can be disrupted by Galactic tidal fields if the internal or- 
bital period is longer than a Galactic year, Pc — 230 Myr. The 
critical semi-major axis is therefore, 

( m P 2 V /3 
a crit = AU %] ~2pc, (Galactic Tides). (9) 

V M ©yr/ 

Since these critical values are not much larger than the 
largest projected separation in our sample, one might be con- 
cerned that the effects of disk and Galactic tides are non- 
negligible. However, the observational data show that the 
binary distribution is well described by a single power law 
within the limits of 3 "5 < A9 < 900", and we are being con- 
servative to ignore both tides because incorporating their ef- 
fects would only serve to place more stringent constraints on 
the perturbers. On the other hand, if the binary distribution 
had shown a power-law break, which would have been indica- 
tive of dark matter, it would have been necessary to investigate 
how much of this signature was actually due to tides. 

4.3. Giant Molecular Clouds 

With typical masses 10 5 — 10 6 M Q and surface density 
£ ~ 5 Mq pc~ 2 , corresponding to a mean density p~6x 
10~ 4 M Q pc" 3 over typical halo orbits, giant molecular clouds 
(GMCs) are clearly in the tidal regime. Because of their low 
density p ~ 0.07pn, they would add very modestly to the ef- 
fect of halo tidal perturbers. The effect of GMCs is further di- 
minished by the fact that /? m ; n (eq. Q) is substantially smaller 
than a typical cloud (and this remains so even if one consid- 
ers the individual "clumps" inside a GMC). Hence, we ignore 
GMCs. As with tides, the effect of doing so is conservative: 
taking account of GMCs would only strengthen any limits we 
obtain. 



4.4. Ionized Binaries 

Even though binary systems are disrupted by perturbers, the 
stars of former binary members escape the system with ex- 
tremely low escape energies due to the interplay of diffusive 
processes and the tidal boundary rather than disappearing im- 
mediately after disruption. As shown by the spectacular im- 
ages of the globular cluster Palomar 5 suffering on-going tidal 
disruptions ( Odenkir chen et alJl2003l) . the mass in a tidal tail 
can be larger than the mass in the parent body. By analogy, 
ionized binaries could in principle stay well inside the obser- 
vational range of interest. 

However, globular cluster escapees are no longer perturbed 
and so proceed on deterministic orbits that separate from the 
cluster at a roughly uniform rate, while the binary members 
drift away in the Galactic potential and continue to be heated 
by the omnipresent perturbers after they are well enough sep- 
arated to be free from each others' gravitational influence. 

In the Appendix, we quantify this argument and analytically 
show that ionized binaries do not contribute significantly to 
the observed distribution. We also give a prescription for in- 
cluding them in our simulations, and indeed we find that they 
contribute negligibly. 

4.5. Non-Circular Orbits 

In our simulations, we assume that the binaries are sub- 
jected to perturbations by ambient objects whose density is 
constant in time. Strictly speaking this would apply only to bi- 
naries in circular Galactic orbits. More typical binaries would 
tend to spend the majority of their time farther from the Galac- 
tic center than the Sun, where the perturber density is lower. 
However, we find by numerically integrating a wide range of 
orbits, that the average density of perturbers encountered by 
binaries in non-circular orbit is almost always higher than for 
those in circular orbits. Hence, our approximation of constant 
perturber density is conservative. 

4.6. Monte Carlo Algorithm 

To test the consistency of the observed binary distribu- 
tion with the presence of massive perturbers, we must be 
able to simulate the effect of these perturbers on arbitrary 
power-law initial distributions. Our basic method for doing 
this is Monte Carlo simulations in which the initial binary 
semi-major axis is drawn uniformly from power-law distri- 
bution over the interval, 1 < log(a/AU) < 5.5. The eccen- 
tricity is drawn randomly from a distribution uniform in e 2 
(Weinberg. Sha piro & Wassermanlll987l) . and the phase and 
orientation of the orbit are assigned randomly. The per- 
turbers are assumed to have an isotropic Maxwellian veloc- 
ity distribution relative to the binary center of mass, with 
a one-dimensional dispersion a = 200 km s" 1 . This repro- 
duces the true rms velocity, which is a combination of an 
isothermal-sphere perturber distribution of circular speed v c = 
\Phi = 220 km s" 1 , and the measured dispersions of halo 
stars 1 (a y, erg, ay) = (170,97,93) km s" 1 ( Popows ki & Gould! 
1998a b). In principle, the velocity distribution should be 
treated as anisotropic. However, this is a higher-order effect, 
which we ignore in the interest of simplicity. 

We consider all impact parameters with b < b max = 
max{10b m in,2a}. In particular, we evaluate the ~ 100 clos- 

1 Chiba & Beers i2000) obtained dispersions of halo stars (a^ ,erg ,a z ) = 
(141,106,95) kms~'. The corresponding one-dimensional dispersion for 
perturbers is a = 194 km s , which is therefore insensitive to the precise 
choice of (<t„- , erg , a-). 
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FIG. 2. — Binary distributions as a function of semi-major axis. 100,000 
binaries are generated following an arbitrarily chosen flat (a=l) distribution 
represented as a thick solid line. The halo density is set to be pn ■ The squares, 
triangles and circles represent binary distributions for three different masses 
of perturber, after T = 10 Gyrs evolution. The fitting curves for each model 
are shown as dashed lines. 



est impacts in the tidal regime, even though as we discuss in 
§13 the single closest encounter dominates. The mass of each 
binary component is set to be 0.5 M Q . We evolve 100,000 
binary systems in each simulation. To illustrate the depen- 
dence on the mass of perturbers, we show their effects on an 
artificial initial binarydistribution that is independent of semi- 
major axis (see Fig.|2j. As expected, for the binary systems 
that are initially tightly bound, the final distributions are al- 
most the same as the initial ones regardless of the mass of the 
perturbers. However, at wide separations, the distributions 
are driven to a new power law, which becomes steeper with 
increasing mass. 

This approach works well for any individual initial power- 
law distribution, but is too time-consuming to process the very 
large number of distrib utio ns required for the comparison of 
data and models. In § 15.51 we introduce a scattering-matrix 
formalism that substantially improves the efficiency of mass- 
production Monte Carlo simulations. 

4.7. Fitting Formula 

Motivated by the fact that the final distribution is well- 
approximated by power laws at each extreme, we use a two- 
line (five-parameter) fitting formula given by, 

H(x)=[f(xr+ g (xr]~ l/n , do 

where f(x) and g(x) are (two-parameter) straight lines in their 
argument, x = log(a / AU), each corresponding to the respec- 
tive asymptotic behaviors of H(x). The fifth parameter n per- 
mits a smooth transition between f(x) and g(x) in the inter- 
mediate region. We calculate the five parameters of equa- 
tion i llOt by minimizing \ 2 f° r a given data set, and the fitting 
curves for four different masses of perturber are represented 
as dashed lines in Figure|2| 



5. RESULTS 

Here, we present our main calculations on the evolution of 
binary distributions with various initial slopes under the influ- 
ence of various perturber masses and halo densities, and we 
evaluate the transition separations. Although the semi-major 
axis of a binary system is a direct indicator of the binding 
energy of the system and is the theoretically most tractable 
quantity, it is not observable. It is the angular separations on 
the sky that we can directly measure from observations. To 
compare our results with the data, we calculate physical sep- 
arations projected on the sky plane and convolve these with 
an adopted distance distribution to predict the binary distribu- 
tions as a function of angular separation. 

In principle, one could compare models directly to the ob- 
served projected separations since CG give individual dis- 
tance estimates to each binary. However, while the observa- 
tional selection function is quite simple for angular separa- 
tions (essentially just a pair of 9-functions), it is rather com- 
plex and would be extremely difficult to model for projected 
physical separations. Hence, we compare our models to the 
most directly observed quantity: angular separations. 

5.1. Semi-Major Axis, a 

We begin the simulation with a flat distribution, 
dN/dloga oc const for calculation of scattering matri- 
ces. We then investigate the dependence of the resulting 
final binary distribution on perturber mass and halo density. 
Figure [3] shows a sample of our results with the initial 
power-law distribution, dN/dloga oc a" 0,567 , in agreement 
with the observations as summarized in §|2] Two models with 
the same perturber mass but widely different halo densities, 
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FIG. 3. — Evolution of binary distributions with the initial power law, a = 
1.567 obtained in §|2| represented as a thick solid line. 50,000 binaries are 
evolved in the presence of perturbers of mass 1000 Mq The circles show the 
binary distributions as functions of semi-major axis while the triangles show 
the distributions of projected physical separations of the binary components 
onto the sky plane. The solid lines represent fitting curves for semi-major 
axis, and the dashed lines for projected physical separation. Two vertical 
arrows indicate transitions, a,, from the unperturbed to the perturbed regimes. 
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Ph and 10 1 pn are shown. The distributions of semi-major 
axis are represented by circles. 

5.2. Projected Physical Separation, r±_ 

At T = 10 Gyr, we calculate the projected physical separa- 
tion, r± of each binary taking account of the randomly chosen 
orientation to the line of sight and the orbital phase. The dis- 
tributions of r± are shown by triangles in Figure|5] 

The two sets of distributions are very similar to each other 
as demonstrated by the solid and the dashed lines, which are 
fitted to the points. This can be understood as follows. The 
time-averaged physical separation is 

/ 2 - 



(r(e)) 



1 

2^ 



a(l — e cos tp) dip = a[ 1 + 



(11) 



Averaging over the uniform distribution in e 2 yields (r) = 
5a / 4, and projecting onto the sky plane gives, 

(r x ) = |(r) = ^a~0.98a. (12) 

That is, the triangles are, on average, slightly shifted to left of 
the circles in Figure|3] 

There is, however, a deviation from equation (11 21 at large 
projected physical separation. We find that a binary system 
with a highly eccentric orbit is more likely to be disturbed 
than one with a circular orbit. Since the highly eccentric sys- 
tems are selectively destroyed for loosely bound binary sys- 
tems, the distribution in eccentricity is no longer isotropic, 
and the actual projected physical separations differ from equa- 
tion dl2> . We therefore use numerical calculations of the pro- 
jected physical separations rather than the analytic estimate of 
equation (I12> . 

5.3. Transition Separation, a t 

Using the five-parameter fitting formula from equation (1101 . 
we calculate the asymptotic slopes for distributions in semi- 
major axis. These are represented as solid curves in Figure[5] 
The intersection of these asymptotic slopes characterizes the 
transition a, from the unperturbed to the perturbed regime. 
The two arrows in Figure|3]represent the intersections for two 
different halo densities. As predicted by equation 0, the tran- 
sition for a denser halo model occurs at smaller separation, 
a, oc jO -2 / 3 . 

The transition separations for various perturber mass and 
halo densities are shown in Figure^] which provides a qualita- 
tive understanding of the underlying physics of binary disrup- 
tion and is in a good accord with the prediction of equation 
and (|6}. In particular, in the tidal regime the transition separa- 
tion is independent of perturber mass, and is a function of halo 
density a, oc p~ 2 ^, although the actual values of a, are a factor 
of 2 smaller than those predicted by equation (|3}. At smaller 
mass, the transition separation grows roughly as a, oc M~ l , as 
predicted by equation for the Coulomb regime. 

5.4. Angular Separation, A9 

To predict the observed angular distribution function, we 
must convolve the model's projected-separation distribution 
function with an assumed distribution of (inverse) distances 
to the binaries in the sample. 

We derive this distribution from the observed distances of 
the 90 binaries in the CG sample. To understand our proce- 
dure for doing so, consider first the idealized case of bina- 
ries drawn from a distance-limited sample of uniform den- 
sity. Of course, in this case, the distance distribution would 
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FIG. 4. — Transition separations a, for different halo densities as a function 
of perturber mass. Transition separations are ob tain ed for each model by 
calculating the intersection of f(x) and g(x) in eq. 1101 . 



be dN/dlogR oc R 3 Q(R mdx - R), where 9 is a step function, 
and R mdx is the distance limit. However, the binaries selected 
at fixed A9 would not be distributed as R 3 . This is because the 
binaries at distances Ri and R2 would have projected separa- 
tions r± \ = R\A8 and r±2 = RiAO, and these differ in fre- 
quency (for fixed log intervals) in the ratio (r±j/r±,2) l ~ a , 
where a = 1.567 is the binary-separation power-law slope. 
Hence, the sample distribution would be dN samp i e /dlogR oc 
R 4 ~ a <d(R m!lx - R). Therefore, to obtain an estimate of the 
underlying distribution of distances that yields the observed 
sample distribution, we simply adopt the observed distribu- 
tion, but weight each binary by /? Q_1 . We estimate R for each 
binary by applying the fourth-order color-magnitude relation 
of CG to the brighter component, except when that component 
has no /-band data, in which case we use the fainter compo- 
nent. 

5.5. Scattering Matrix, S 

To simultaneously investigate large subsets of initial power- 
law distributions for each model s pecifi ed by M and p, we 
use a Monte Carlo simulation (see § 14. 6i with an initially flat 
(dN /dloga oc const) distribution to construct a scattering ma- 
trix, Sy(M, p), the relative probability that a binary with initial 
log semi-major axis, loga ; will finally be observed in the an- 
gular separation bin, log Af?,. Explicitly, 

S (M,p) = ^ I dlog/?-T^<f[log*A0 / -logr ±J (*;M,p)], 

k °^ 

(13) 

where r± j(k;M, p) is the final projected separation of a binary 
with initial semi-major axis aj in k ,h realization of a Monte 
Carlo simulation of perturbers with mass M and density p. In 
practice, since the radial profile dN/dlogR is e stimated di- 
rectly from the binary data as discussed in § 15.41 we evaluate 
equation Jl 3I > by discretely summing over the / = 1 . . . 90 bi- 
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FIG. 5. — The best-fit final binary distributions for various perturber masses, 
assuming that the initial distribution is a power-law. The halo density is set 
to be pn- The observed halo binary distribution (Fig. is shown for com- 
parison. A model with 1000 Mq perturber deviates significantly from the 
observations while a model with 10 Mq is quite consistent with them. 
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|log^A(9,-logr ±J ^;M,p)| 



where <51og AO is the logarithmic width of the AO bins. 

Since the simulations are calculated from a flat distribution, 
i.e., equal number of binaries per logarithmic semi-major axis 
bin, the un-normalized probability of finding binaries at AO 
for a given dark-matter model is then given by, 



P( A6>,;M, p, a) = S U (M, p)w{oj\ a), 



(15) 



where w(aj;a) = a l j a is the power-law of the initial binary 
distribution, dN /da oc a~ a . 

6. DARK MATTER LIMITS 

For each model, we fix the normalization so that the number 
expected in the interval 3. "5 < AO < 900" is equal to 90, i.e., 
the total number observed over this interval. We call the nor- 
malized resulting function P^(A6;M, p,a), evaluate the like- 
lihood of a given model as a sum over the 90 observed bina- 
ries, 

90 

In L(M, p,a) = J2 lnP N (A6»,;M, p, a), (16) 
;=i 

and find the maximum likelihood over various initial slopes 

a, 

In L( M, p) = max {In L(M,p, a)}. (17) 



Figure |3 shows the normalized resulting functions for var- 
ious perturber masses with p = pn- For 1000 M Q perturbers, 
the likelihood is maximized by an initial distribution that is 
nearly fiat. However, no initial power-law is consistent with 
the observed distribution. For low mass perturbers, by con- 
trast, the likelihood is maximized by an initial slope that is 
only slightly shallower than the observed one, and the result- 
ing distribution is consistent with observations. 

We compare the resulting likelihoods to L(0,0), i.e., the 
model with no dark-matter perturbers, which is characterized 
by the best fit slope a = 1 .567, and we define confidence levels 

by 

a(M,p)= V2[lnL(0,0)-lnL(M,p)]. (18) 

As we discussed in § 12 the fact that the binary sample 
is clearly incomplete for AO < 3 "5 implies that it may also 
be slightly incomplete for binaries just above this threshold. 
However, one can see from the form of the normalized re- 
sulting functions in Figure[5]that such incompleteness would 
actually favor models with massive perturbers over the no- 
perturber model. We verify numerically that this indeed is the 
sign of the bias. Hence our choice of the AO > 3 "5 threshold 
is conservative in the face of possible incompleteness. 

Figure [6] shows the resulting contour plot excluding vari- 
ous dark-matter models at various confidence (a) levels. Fi- 
nally, in Figure0we comp are our 95% confide nce (2a) limits 
with those from the EROS (lAfonso et alJ2003l) and MACHO 
jAlcocketalJl200l microlensing collaborations. The halo 
binaries exclude models that have generally higher MACHO 
masses than those probed by these microle nsing ex periments. 
Moreo ver, they extend all the way to the lLacev & Ostrikerl 
(1985) limit based on stability of the Galactic disk. 

All three curves in Figure0show limits on models with 5- 
function mass distributions. However any model with broader 
distribution of mass but total density pu, would be ruled out 
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FIG. 6. — Exclusion contour plot for halo dark-matter models. Various con- 
fidence levels are shown. The oscillations at high masses are due to numerical 
noise. 
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FIG. 7. — Exclusion contour plot at 95% confidence level. The dashed, 
the dotted and the long dashed lines represent the microlensing-based limits 
from EROS i Afonso et al. 2003;) and MACHO lAlcock et all2001l) . and the 
limit based on disk stability ILac ev & Ostrikei 1 985j), respectively. Our limits 
from the 3" 5 < A9 < 900" sample are represented as a solid line. 



7. CONCLUSIONS 

In this paper, we have investigated the evolution of halo 
wide binaries in the presence of MACHOs and estimated 
upper limits of MACHO density as a function of their as- 
sumed mass by comparing our simulations to the sample of 
halo wide binaries of CG. We exclude MACHOs with masses 
M > 43 M Q at the standard local halo density p H at the 95% 
confidence level. 

MACHOs have been a major dark-matter candidate ever 
since observations first established that this mysterious sub- 
stance dominates the mass of galaxies. Prodigious efforts over 
several decades have gradually whittled down the mass range 
allowed to this dark-matter candidate. However, the window 
for MACHOs with 30 M Q < M < 10 3 M Q remained com- 
pletely open while constraints in the range 10 3 M < M < 
10 6 M Q were somewhat model dependent. Our new limits on 
MACHOs M > 43 M all but close this window. 



provided that the combined limits from the three curves ruled 
out each mass separately. 

Hence, the only remaining window open for MACHOs 
would be black holes with a mass function strongly peaked 
atM ~ 35 M Q . Even this model is more strongly constrained 
than shown Figure because the MACHO lAlcock et alJ 
120011) results and our results each separately place weak limits 
on this model, which when combined comes close to ruling it 
out. 
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APPENDIX 
IONIZED BINARIES 

After disruption of a binary system, the two ionized members remain in similar Galactic orbit, so it is only the separation along 
the direction of the orbital motion that can keep increasing, while the perpendicular separation oscillates. 

In the Coulomb regime, the average post-ionization gain in the relative velocity of binaries along the orbital direction is (see 
eq.0), 

where At is the remaining time to 10 Gyr after disruption. For a diffusive process that is uniform over time At, the root-mean- 
square separation parallel to the orbital motion is 

where the Coulomb logarithm is calculated at the scaled quantities. 

For a conservative limit of the tidal radius, a, = 3 pc, the time required for ionized binaries to separate farther than a, is 0. 1 3 Gyr, 
so that only binaries ionizing within the last ~ 1 % of the age of the Galaxy have a significant chance to be confused with bound 
systems, even for the lowest mass perturbers that we can effectively probe. In the tidal limit, binaries escape with characteristic 
velocities of the transition separation ~ 10 4 AU, i.e., 300 m s" 1 . Hence, they drift one tidal radius in only 10 Myr, so their impact 
is even smaller. Nevertheless, since there are only a handful of binaries in the widest-separation bins, it is important to make a 
careful estimate of the contribution from ionizing binaries. We take account of ionized binaries in the simulations as follows. 
Binaries are considered ionized either when they have positive energy or they have a > a,. They are then assigned a relative 
velocity equal to their escape velocity in the former case, or zero in the latter. A random orbital direction is chosen. The ionized 
binaries in the Cou lomb regi me then continue to suffer perturbations to the end of the simulations whose effect we calculate 
using equation (IA1> and (IA2> . At the end of the simulation, the binary is assigned a transverse separation drawn randomly from 
a sinusoidal distribution of amplitude, 

^MAX=^=26pc(^i rr ), (A3) 

where vj_ is the final transverse velocity, Q = ^/2v c /Rq is the epicyclic frequency, and Rq is the Galactocentric distance. Finally, 
we "observe" the ionized binary from a random orientation and record the projected separation. We find that the ionized binaries 
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have no significant effect either on the final distribution or the calculated likelihoods. We ignore ionized binaries in the tidal 
regime because these esca pe with much higher initial velocities and because these velocities grow much more rapidly than 
indicated by equation lAlll . 
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